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Abstract 

We report the results of computer simulations of epitaxial growth in the presence of a 
large Schwoebel barrier on different crystal surfaces: simple cubic(OOl), bcc(OOl), simple 
hexagonal(OOl) and hcp(OOl). We find, that mounds coarse by a step edge diffusion driven 
process, if adatoms can diffuse relatively far along step edges without being hindered by 
kink-edge diffusion barriers. This yields the scaling exponents a = 1, j3 = 1/3. These 
exponents are independent of the symmetry of the crystal surface. The crystal lattice, 
however, has strong effects on the morphology of the mounds, which are by no means 
restricted to trivial symmetry effects: while we observe pyramidal shapes on the simple 
lattices, on bcc and hep there are two fundamentally different classes of mounds, which are 
encompanied by characteristic diffusion currents: a metastable one with rounded corners, 
and an actively coarsening configuration, which breaks the symmetry given by the crystal 
surface. 

1 Introduction 

In spite of considerable efforts [1-7], see || for an overview, a thorough theoretical under- 
standing of the late phases of epitaxial crystal growth is still lacking. In this publication we 
investigate the problem of growth in the presence of a strong Schwoebel barrier, which hinders 
interlayer transport and leads to an instability of the flat surface. During an initial transient, 
mounds form on the surface, which then start to merge. It is generally accepted, that in this 
asymptotic coarsening regime the statistical properties of the surface remain invariant under 
a simultaneous transformation of spatial extension x, height h(x) and time t: 

x -» bx; h(x) -» b a h; t -» b z t, (1) 
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where a and (3 := a/z are believed to be universal exponents, which do not depend on details 
of the model and b is an arbitrary factor. If the process of growth stabilizes a specific slope, 
this will yield a = 1. 

It was first pointed out by Siegert et. al. [||, ||], that lattice symmetries may play an 
important role in the coarsening process. These authors investigated continuum equations, 
using an analogy between coarsening and a phase ordering process, which has recently gained 
popularity ||: areas of constant slope should correspond to domains of a constant order 
parameter. They derived scaling exponents a = l,/3 = 1/3 on surfaces with a triangular 
symmetry, and j3 = l/(3\/2) for generic cubic surfaces, while (3 = 1/3 requires a fine-tuning 
of parameters. 

However, Monte-Carlo simulations [10, |ll|] have raised doubts on these predictions, since 
they yield (3 ~ 1/3 on a simple cubic lattice for a great range of parameters. In this paper, we 
want to adress the following questions: (1) What are the mesoscopic processes which make 
the mounds coarse? (2) How does the coarsening process depend on the crystal lattice and 
its symmetries? (3) Do these results support a deeper analogy between coarsening and phase 
ordering? 



2 The model 

To answer these questions, we perform computer simulations of growth on (OOl)-surfaces of 
the simple cubic (sc) lattice, the simple hexagonal lattice (sh), the body-centered cubic (bcc) 
lattice and the hexagonal close packing (hep). This is done under solid-on-solid conditions, 
i.e. the effects of overhangs or dislocations are being neglected. Then, the simple lattices 
can be represented by a square (sc) respectively a triangular mesh (sh) of integers, which 
denote the height h(x) of the surface. We build the bcc (hep) lattice out of two intersecting 
sc (sh) sublattices, which contain the even respectively odd heights^]. Here, an adatom in a 
stable configuration is bound to 4 (3) neighbours below in the other sublattice, which will be 
denoted as vertical neighbours in the following. Particles with fewer vertical neighbours form 
overhangs, which are forbidden by the solid-on-solid condition. This is physically reasonable, 
since such particles are only weakly bound and therefore these configurations will be unstable. 
Additionally it may have neighbours in the lateral direction, which are in the same sublattice. 

The investigation of the coarsening process requires a fast algorithm which allows for 
the simulation of the deposition of thick films on comparatively large systems. Standard 
kinetical Monte-Carlo techniques, which consider the moves of many particles on the surface 
simultaneously require computationally expensive bookkeeping procedures. This makes them 
too slow for our purpose. Instead, we simulate the moves of a single particle from deposition 
until an immobile state is reached |u], Q| : 

An adatom impinges on a randomly chosen lattice site. Due to its momentum perpen- 



dicular to the surface, the particle may funnel downhill [12, 13, 14] to the lowest (vertical) 
neighbour site. On bcc and hep this is repeated until it reaches a stable site, as defined above. 
On the simple lattices this is a site x, where all nearest neighbour sites have a height > h(x). 

Then, the adatom diffuses on the surface. If the particle has no lateral neighbours, one 
of its neighbour sites is chosen at random. On the simple lattices, the particle moves to this 
site only if its height does not change, i.e. we introduce an infinite Schwoebel barrier. On the 



1 For simplicity, we assume the spacing between the layers to be 1 lattice constant. Our algorithm depends 
only on the topology of the lattice. 
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bcc or hep lattice, the particle is moved to the neighbour site if it is stable. This condition 
implies an infinite Schwoebel barrier, too, since a transition into a different layer can proceed 
only via a weakly bound (unstable) transient state. This procedure is repeated until either a 
lateral neighbour is reached or 1^ steps have been performed. The simulation of single particles 
requires an effective representation of the collisions of diffusing adatoms, which form a stable 
nucleaus. The diffusion length corresponds to the mean free path length and depends on 
the diffusion constant D and the particle flux F: oc (D/F) 1 ^ |l5|| . 

As soon as a particle has a lateral neighbour, it is bound to it. This process is irreversible 
in the sense that we forbid diffusion processes which reduce the number of bonds, like the 
detachement from a step edge. However, the adatom may diffuse along the edge. After 
l\ steps, or if the particle has reached a kink site, it is fixed to the surface. If not stated 
otherwise, diffusion around corners is allowed. 

Within this model, we measure time t in units of the time needed to deposit a monolayer. 
This algorithm may be programmed very efficiently and is about an order of magnitude faster 
than full diffusion Monte-Carlo algorithms. 

In all our simulations, we choose Id = 15. The diffusion length sets the initial island 
distance only, while its influence can be neglected in the later phases of growth, when the 
typical terrace width is much smaller than l^- We simulate a variety of different values of the 
step edge diffusion length l^ in the range between 1^ = 1 and 1^ = 20. The simulations are 
performed on a square oi N x N lattice constants using periodic boundary conditions (bcc 
and sc), respectively a regular hexagon with edges of length M (hep and sh) using helical 
boundary conditions, our standard values being N = 512 and M = 300. For every parameter 
set, we have performed 7 independent simulation runs. 



3 Surface morphologies 

During the first ~ 100 monolayers of growth on an initially flat surface, islands nucleate, on 
which mounds build up and take on their stable slope. Then, the asymptotical coarsening 
regime starts. As expected, on the simple lattices the mounds obtain regular shapes, which 
are determined by the symmetry of the surface: square pyramids on sc, hexagonal ones on 
sh. On bcc and hep however, rounded corners as well as sharp corners can be found (figures 
|] b, §). Here one finds two fundamentally different types of mounds: On the one hand, 
there are mounds, where all corners are rounded with octagonal shapes on bcc (figure 2a), 
and 12-cornered ones on hep. On the other hand, one observes a breaking of the symmetry 
given by the lattice: approximately triangular mounds on hep (figure || b) and oval shapes 
on bcc (figure ||b), where every second of the corners is sharp, while the others are rounded. 
We have observed, that in the process of coarsening, on all lattice structures mounds merge 
preferentially at the corners, while mounds touching at the flanks are a metastable configu- 
ration. On bcc and hep, the metastable configuration is made of mounds with completely 
rounded corners. Mounds, whose shapes break the lattice symmetry, are always in the course 



of merging with a neighbour they touch with a sharp corner [16]. 



To characterize the morphology quantitatively, we investigate the slopes m = V/i of the 



surface. On average, |m| obtains a value which can be estimated by a simple argument [17]: 
The Schwoebel effect yields an uphill current due to the preferential attachement of particles 
on terraces from below, which is compensated by a downhill current due to funneling. If we 
assume, that a particle reaches its final height after at most one incorporation step, the sc and 
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Figure 1: Simulation results on hep (M = 300, l d = 15, l k = 10, t = 20000ML). Panel a: 
Histogram of slopes (density plot). High probabilities are drawn dark. Panel b: Contour plot 
of a part of size 250 x 250 lattice constants of the surface. High levels are plotted in light grey. 
The arrows show the surface current, as measured during deposition of additional 200ML. 

the sh surface select |m| = 1/2, while bcc and hep select \rh\ = 1. In practice, the selected 
slopes are slightly smaller, since the maximal number of incorporation steps in the funneling 
process is unlimited. Since a direct computation of the numerical gradient of simulated 
surfaces fails due to the discrete heights, we first apply a gaussian filter with variance a = 4. 
This value has turned out to be a good compromise: it removes the atomistic structure, 
but preserves the shape of the mounds. Figure |l|a shows a two-dimensional histogram of 
the slopes on hep surfaces: due to the round shapes of the mounds, one finds a pronounced 
maximum of the probability density function in every direction of fh. The most probable value 
of \rh\ has a directional dependence: it is minimal in the lattice directions (0.88 ± 0.02), and 
maximal (0.99 ± 0.02) in the intermediary directions. Consequently, the rounded corners are 
steeper than the flanks. On bcc the results are equivalent, apart from the different symmetry, 
while they are rather trivial on sc and sh: a regular square respectively hexagon of slightly 
broadened peaks, which correspond to the flanks of the mounds. 

4 Diffusion currents 

To gain deeper insight into the coarsening process, we investigate the material transport on 
mesoscopic lengthscales, which can be done by tracing the motion of particles on the surface. 
On all lattice structures, we observe a strong uphill current at the corners of mounds, which 
is a geometrical effect: an adatom attaching at the corner is transported over a mean distance 
lk, namely with equal probability along either of the two flanks of the mound. This results 
in a net current j in the direction of the bisector of the angle at the corner, which is directed 
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Figure 2: Contour plots of parts of size 250 x 250 lattice constants of bcc surfaces (N = 512, 
Id = 15, Ik = 10, t = 20000ML). Panel a: section dominated by quickly coarsening mounds 
breaking the surface symmetry. Panel b: metastable configuration. 

uphill. A simple calculation yields 

\j\ = r a l k cos - 

if there is a straight step edge of length S> Ik, which is the case on the simple lattices. Here 
r a is the rate paricles attache at a step edge with, (ft = it/2 on surfaces with cubic symmetry 
and (j) = 27r/3 on hexagonal surfaces. 

On the bcc and the hep lattice however, material transport along step edges is severely 
restricted due to the rounded corners. There, particles diffusing along the step edge are on 
average reflected, which leads to a current towards the middle of the smooth step edges at 
the flanks of mounds in the metastable configuration, which have an extension ~ Ik (figure 
2b). This current has a non- vanishing uphill component due to the Schwoebel effect. On the 
rounded corners on the contrary, one observes a weak downhill current. This is explained by 
the observation of steep gradients (figure [l] a), at which the downhill current due to funneling 
dominates the uphill current of the Schwoebel effect. On these surfaces the stable slope, 
at which the currents cancel, is not assumed locally, but only in the global average, which 
results in spatial current patterns on the surface. An analogous behaviour is observed at 
the symmetry-breaking, merging mounds (figures lb, 2a). Here, the smooth edges are in the 
proximity of the sharp corners. Particles are reflected at the rounded corners, too. This yields 
a net current towards the sharp corners, which compensates the geometrically induced uphill 
current exactly, if the smooth edges have a length Ik- This fact explains the broken symmetry 
of these mounds: if there is any process, which transports matter towards sharp corners, this 
material is not completely diffusing inward, as would be the case on pyramidal mounds, but 
is attached to the corner and makes it overgrow its neighbours. Merging of mounds, however, 



is such a process: consider a terrace, which goes around two mounds touching each other at 
the corners. Due to its curvature, in this contact zone there is a high concentration of kink 
sites, which make it a sink for diffusing particles. 



5 Theoretical estimation of scaling exponents 

Our observation, that mounds merge perferentially in positions with touching corners, gives 
strong evidence that the current which fills the gap between two mounds, plays a dominant 
role in the process of coarsening, at least in the case of large l^, where it is efficient. If 
this is the case, then a simple consideration yields the scaling exponents: Due to step edge 
diffusion, adatoms are transported over a typical distance l^. If 1^ is noticeably greater than 
a few lattice constants, this is much larger than the average terrace width. Then, material 
transport via diffusion on terraces is small compared to transport via step edge diffusion and 
can be neglected. Figure |l]b shows, that the current into the gap is significant on a few atomic 
layers in the vicinity of the contact point only. In consequence, this current is independent of 
the size of the mounds, if the latter is much greater then l^. Since the volume of the gap is 
proportional to L 3 , if L is the typical distance between the mounds, it will take a time t ~ L 3 
to fill it. In consequence, L ~ i 1 / 3 ==>• z = 3. Since a = 1 in the presence of slope selection, 
this yields (3 = 1/3. In contrast to the effects discussed in ||, ||] these considerations are 
completely independent of lattice symmetries. In Ij, g a similar argument has been applied 
to the case of coarsening in the absence of dominant step edge diffusion processes ("bond 
energy driven coarsening"), where material is transported mainly by terrace diffusion. This 
yields j3 = 1/4. The same exponent is obtained, if coarsening is exclusively due to fluctuations 
in the particle beam ("noise assisted coarsening"). We expect, that these processes dominate 
step edge diffusion in the limit of small 4, which leads to a transient from /3 = l/4to/3 = l/3, 
when Ik is increased. 

6 Measurement of scaling exponents 

We apply a variety of methods to measure the scaling exponents. This is important as a consis- 
tency check and to eliminate systematical errors which might be intrinsical to some methods. 

/ - 2V 1 / 2 

First, can be obtained from the increase of the surface with w(t) = ( (h(x, t) — (h) (t)) ) 

with time. Equation [l] corresponds to w(t) ~ t@ ]l8|]. This power law behaviour may be cor- 
rupted by the presence of noise on the surface profile, an aditional contribution to w(t) which 
has been called intrinsic width Similarly, the number of mounds n m in the system will 
decrease like t~ 2 ^ z , if this scaling hypothesis holds. We measure n m by counting the number 
of top terraces in the system. Since a single particle is counted as a terrace, this method 
may be misleading, if particles nucleate on terraces at the flanks of mounds. A method 
which is widely used in the literature uses the fourier transform /i(fc)0: The structure factor 
S(k) := (h(k)h(—k)^j will be maximum at nonzero wavenumbers \k m \ = 2n/l m , if there are 

structures at a typical lengthscale l m on the surface. Since l m ~ t 1 ^, \k m \ will decay ~ t~ x l z . 
In practice, a direct search for the maximum often fails due to noise effects; one avoids this 
problem by calculating the averages km^ •= ((j2k\k\S(k) p ^j /Ej:S(k) p ), however, the choice 

2 The mean surface height is subtracted before performing the fourier transform. 
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of the correct power p is a bit arbitrary. 



6.1 The wavelet method 

To avoid these difficulties, we propose a method which exploits the continuous wavelet trans- 
form 

f^]M:=i/^f ^fi^ W) (2) 
a J \ &2{(x - b)/a) J 

of the surface heights. Here, the wavelets \&2 are defined as partial derivatives of a radially 
symmetrical filter which we choose to be a gaussian: ^1 := d^/dx, ^2 '■= d<&/dy. The 

basic idea ot this transform is to convolve h(x) with the wavelets, which are dilated with the 
scale a. This makes it a natural tool to search for the typical lengthscale at of structures on the 
surface. In |2(], |2l|], it has been proposed to investigate the wavelet transform modulus maxima 
(WTMM). They are defined as local maxima of the modulus Mx[h](b,a) := |T^[/i](6, a)\ in 
the direction of T^. The search for the WTMM is equivalent to edge detection: they lie on 
connected curves, which trace the contours of objects of size ~ a. We investigate the average 
of the WTMM on the scale a 

W m (a):=(M^[h](b,a)) i;=WTMM . (3) 

This function has a pronounced maximum at a value a m tx at, which is the scale that contains 
the most relevant contributions to the surface morphology. On mounded surfaces, a m ~ t 1 ^ 
is proportional to the lateral size of the mounds, while W m (a is a measure for their 

heights. This procedure has several advantages compared to the standard methods used in 
the literature: (1) The detection of the WTMM leads to an efficient suppression of noise, 
and therefore eliminates the intrinsic width. (2) Since only the most important lengthscale is 
considered, the results may not be corrupted by details of the surface morphology on small 
lengthscales. (3) Appropriate wavelets are well localized both in real space and in fourier 
space. Therefore, this analysis combines the advantages of both real space techniques, like 
the counting on mounds, and fourier techniques, like the calculation of k^. We have tested it 
both with artifical test surfaces and various toy models of growth and found a clear superiority 
to the standard methods, especially in the presence of noise, which should make it interesting 
for experimental investigations. Details will be given in a forthcoming publication. 



6.2 Results 

In our simulations, the asymptotic scaling behaviour is obtained after a long initial transient 
of ~ 100 monolayers. Since after this transient only a comparatively small number of mounds 
is left on the surface, the measurement of exponents is complicated by finite size effects: due 
to the periodicity of the system, there is a self-interaction of the currents on the mounds, if 
their size is no more small compared to the system size. We measure a ~ 1 for all Ik > 1, 
a consequence of slope selection. At Ik = 1 however, we measure smaller values, which do 
depend on the crystal lattice: a = 0.85 on bcc and a = 0.94 on hep. A similar effect was 
observed in |lC| , |TTj| on the simple cubic lattice. We remark, that our results are independent 
of the method which was applied to measure the exponents. In figure 3, (3 is plotted as 
a function of Ik- Clearly, its value does not depend systematically on the symmetry of the 
crystal surface. We find no significant deviations between (3 on bcc and on hep (sc and sh). 
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Figure 3: (3, as obtained from the 
wavelet method, as a function of the 
step edge diffusion length. Diamonds: 
bcc lattice, squares: sc, stars: hep, 
circles: sh. Results from measure- 
ments of (3 from the surface witdh are 
identical within the error bars. These 
have been estimated from the loga- 
rithmic fits; we expect the true errors 
to be larger due to systematical devi- 
ations. 



As already observed in |T(| p] , there is a dependence of (3 on the step edge diffusion length. 
One obtains values « 1/4 at Ik = 1 and higher values at greater 1^. The values obtained 
at large l). are compatible with (3 = 1/3. As indicated above, we explain this transient as 
a competition between different coarsening mechanisms: At large Ik the merging of mounds 
should proceed mainly by the step edge diffusion driven transport of material into the gap, at 
small Ik it is dominated by noise assisted coarsening and/or bond energy driven coarsening. 
This competition might also explain, why the exponents measured on bcc and hep are still 
slightly smaller than those obtained on the simple lattices even at values of Ik as large as 20, 
since the symmetry breaking of the mounds on these lattices halves the number of corners 
which are active in the coarsening process. 

All the simulation results reported above have been obtained in simulations with unhin- 
dered step edge diffusion. To further corroborate our ideas of the coarsening process, we have 
also performed simulations with a corner diffusion barrier. A particle may diffuse at most Ik 
steps along a straight step edge, but diffusion around corners is forbidden. Then, we observe 
no symmetry breaking of the mounds, which now obtain regular polygonal shapes on all lat- 
tice structures. These are rotated by an angle of 45 degrees against the lattice directions on 
cubic surfaces and 30 degrees on hexagonal ones. Due to the absence of aligned step edges, 
material transport by step edge diffusion is severely restricted under these conditions. Thus, 
we measure only small diffusion currents, and observe no pronounced long-range order in the 
flux lines. Here, we find a « 1, /3 ~ 0.25 on all lattices, independent of the step edge diffusion 
length. This slow coarsening in absence of the characteristic features of step edge diffusion 
driven coarsening - currents on mesoscopic lengthscales and symmetry breaking on bcc and 
hep - strongly supports the important role of this mechanism for fast coarsening with (3 = 1/3. 

7 Conclusions 

In summary, we have presented a detailed investigation of the coarsening process in epitaxial 
growth. Our most important finding is, that the crystal lattice has a considerable influence on 
the morphology of the growing surface, which is by no means restricted to trivial symmetry 
effects. In spite of these differences, in the limit of large step edge diffusion lengths, the 
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mounds coarse according to power laws with universal exponents. We obtain (3 ~ 1/3 in 
the case of unhindered step edge diffusion and as 1/4 in the case of restricted step edge 
diffusion. 

These results contradict previous studies of the coarsening process using continuum equa- 
tions Q [|, ||], which predict slower coarsening on surfaces with a cubic symmetry. In these 
publications, equations have been studied, which are invariant under the transformation 
h(x) — > —h{x). This symmetry reflects itself also in the up-down-symmetry of the solu- 
tions of these equations. Clearly, our simulation results break this symmetry. In figures [l|a, || 
the contours of the mounds and those of the valleys are clearly distinct. This is also true for 
the currents, especially the step edge diffusion current, which dominates the coarsening be- 
haviour in the limit of large Ik and unhindered step edge diffusion and determines the scaling 
exponents. Since the symmetrical equations do not only fail to describe the morphology of 
the surface correctly, but also predict the wrong universality classes, we conclude, that the 
breaking of the up-down symmetry is a central feature of unstable epitaxial growth. 

However, we find it difficult to interpret this result in the context of an analogy to a phase 
ordering process, where the local slope of the surface corresponds to an order parameter. In 
this picture, the importance of the breaking of up-down symmetry leads to the conclusion, 
that now the stability of a domain wall between areas of constant slope is a complicated 
function of both the orientation of the wall and the order paramters on both sides of it. 
Additionally, in contrast to the simple lattices, where there are only 4 (sc) respectively 6 (sh) 
stable slopes, on bcc and hep one would have to deal with an order parameter, which obtains 
continuous values. 

In any case, the behaviour of our models, which implement the microscopic processes on 
the surface, is governed by rules, which are much more complex than those implemented in 
all the continuum models we know. Differences are by no means restricted to some details, 
but have a fundamental influence on the behaviour of the system on large lengthscales. 
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